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Abstract 

We propose a model for a two dimensional, associative water-like lattice gas 
with one single variable representing both long and short-range interactions. 
The corresponding hamiltonian was solved exactly, by state enumeration in 
a finite lattice, so to obtain an analytic expression for the partition function. 
The lattice dimensions were chosen based on geometric characteristics of the 
stable phases found in previous works using Monte Carlo simulations. An 
expression for the particle density in the finite lattice was then obtained, 
and coexistence curves with a region of anomaly in the density presented 
in a phase diagram. In the end, a phenomenological theory for the system 
density is proposed and compared to the previous results. 
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1. Introduction 

Although ubiquitous in nature and subject of many studies, water is a 
substance that is not yet completely understood by scientists 0,1, 3- While 
most liquids contract upon cooling, the specific volume of water starts to 
increase at ambient pressure when cooled below T = 4°C 0, 0|. Besides, 
within a certain range of pressures, water also exhibits an anomalous in- 
crease of compressibility and specific heat upon cooling jsl, 0, 0] • One of the 
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hypothesis for the presence of the anomalies is that they are related to a 
second critical point between two liquid phases, a low density liquid (LDL) 
and a high density liquid (HDL) [ij]. 

Water is not an isolated case. There are also other examples of tetrahe- 
drally bonded molecular liquids such as phosphorus 10] 11] and amorphous 
silica 12[ that also are good candidates for having two liquid phases. More- 
over, other materials such as liquid metals 13]| and graphite lij also exhibit 
thermodynamic anomalies. Unfortunately, a coherent and general interpreta- 
tion of the low density liquid and high density liquid phases is still missing. 
What type of potential would be appropriated for describing the tetrahe- 
drally bonded molecular liquids? In order to address the question of the 
potential shape appropriated to describe anomalous systems, a few years ago 
it has been shown that the presence of two liquid phases can be associated 
with a potential having an attractive part and two characteristic short-range 
repulsive distances. The smallest of these two distances is associated with 
the hard core of the molecule, while the largest one with the soft core. In the 
case of lattice models, the main strategy has been to associate the van der 
Waals attraction with the occupational lattice gas variable and the hydrogen 



bond disorder with bond 15, 16 



17, 181 Potts states. In the former 



or site 

case, coexistence between two liquid phases may follow from the presence 
of an order-disorder transition, and a density anomaly is introduced ad hoc 
by the addition to the free energy of a volume term proportional to a Potts 
order parameter. In the second case, it may arise from the competition be- 
tween occupational and Potts variables introduced through a dependency of 
bond strength on local density states. In both cases, however, the connection 
between the number of Potts states and the presence of the density anomaly 
is not yet understood. 

To circumvent this difficulty, recently another type of the associated lat- 
tice gas (ALG) model was proposed EEH- Th e main idea is to represent 
hydrogen bonds through ice variables]!]] H > 

so successful in the de- 
scription of ice jiif entropy, for dense systems. In this case, an order- disorder 
transition is absent. Competition between the filling up of the lattice and 
the formation of an open four-bonded orientational structure is naturally in- 
troduced in terms of the ice bonding variables and no ad hoc introduction 
of density or bond strength variations is needed. This approach bares some 
resemblance to that of some continuous models 26] [27] [28J, which, however, 
lack entropy related to hydrogen distribution on bonds. Also, the reduction 
of phase-space imposed by the lattice allows construction of the full phase 
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diagram, not always possible for continuous models. Henriques et al 19] have 
found that this model exhibits two liquid phases and a density anomalous 
region in the pressure-temperature phase diagram. 

The remainder of the paper goes as follows. In sec.©, a simplified model 
is presented and the simplifications regarding its polarity and degrees of 
freedom explained; the exact finite lattice analysis of the model is shown on 
sec.©; simpler expressions for the density and pressure of the model are 
proposed in sec.fll]), and our findings summarized in sec.©. 



2. The Model 

Many existing models of water-like fluids use polar representations for the 
hydrogen bonds between molecules [20] . As a starting point for our approach, 
we take the hamiltonian of equation [Tj for a system that allows 18 possible 
orientational states per particle and presents a phase diagram with liquid 
and gas phases as well as density anomaly. 
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Where r\ denotes the /-arm of the site i that points to the site j, and 
t™ 1 ' 3 ' 1 is the opposite arm from the site j pointing to arm I of the site i, 
selected by Tn^i, so to guarantee that both arms are allowed to participate 
in the effective bonding. 

Such polar formulation for bond formation requires several evaluations 
to determine if a bond between two molecules has occurred or not. In that 
model, the 18 molecular states generate n 18 + 1 terms in the evaluation of 
the partition function in a lattice with n sites, where the "+1" term is due 
to the degenerated empty site states. 

Because of the complexity of all relevant arrangements among particles 
in models involving long and short-range interactions, mean field approxi- 
mations approaches (e.g. Bethe-Peierls) may eventually fail to capture some 
of the system phases, and most results in the literature are obtained from 



Monte Carlo or Molecular Dynamics simulations [20 



In this work, we propose and discuss a model similar to equation [H mod- 
ified so to allow analytical calculations by exact state enumeration in an 
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approach similar to that of Wu et al 29] and Kawabata et al 3o|. Analyz- 
ing the microscopic structure of the outcoming results from previous Monte 
Carlo simulations, we learned about the geometries of the stable phases and, 
based on these geometries, we set the dimensions for the smallest lattice that 
could accommodate them. Then, an exact summation was done, by enu- 
merating all possible states on the finite lattice, obtaining an approximation 
for the grand-canonical partition function of the system, and validating the 
proposed model. 

According to the term Oi<jjTiTj(\ — TiTj) in equation [TJ for all the possible 
combinations of two neighbor arms, resulting bonds are 



{+-, -+} ->■ bond 

{+ + , - - , 00, 0+, - , - , + 0} ->■ non - bond ' 1 ' 

so that probability P(i,j) for a bond occurring between two candidate 
leg states in adjacent sites i and j is 

P{bond) = 2/9 

and (3) 
P(non — bond) = 1 — P(bond) = 7/9. 

It is important to notice that these are not physical Boltzmann probabil- 
ities, so we are ignoring the corresponding energy costs for now. 

To reduce the number of microstates and facilitate the evaluations of the 
grand partition function, we also did two additional simplifications. The first 
one has to do with considering the effectiveness of the bonding directly in 
the hamiltonian, instead of each molecular "leg" polarity. 

If we change the states for each r "leg" from {—1,0,1} to just {0,1}, 
the system hamiltonian can be rewritten, changing the long-range term from 
(TiCTjTiTjft — T^j) to just — 2(7j(7 J TjTj. The original model geometry, consisting 
of four legs allowed to form bonds and two empty links in the triangular 
lattice was kept unmodified, producing the new set of only 3 states, which 
are basically |7r (120°) rotations of the original X-shaped structure that 
represent molecules on the triangular lattice, as can be seen in figure [TJ 

The new probabilities for a bonding to occur, P', still ignoring the corre- 
sponding energetic costs, are 
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Figure 1: Reduced set of states for 2d liquid on the triangular lattice. Lattice sites 
containing molecule extrema are represented by a number 1, otherwise by a 0. 



P' 



(bond) = P(l, 1) = (4/12)(4/12) = (16/144) = 1/9 



and 



(4) 



P'(non - bond) = 1 - P'(bond) = P(l, 0) + P(0, 1) + P(0, 0) = 8/9. 

As it can be seen from equation HJ the new probabilities of randomly 
generating two connected adjacent bonds, P' = 1/9, compared to the original 
P = 2/9 (equation |3]) , still is much smaller than the probability of not having 
any bonding, allowing the expectation of a behavior similar to that observed 
in the original polar model. 

To conclude our set of proposed modifications, we will show that it is also 
possible to get rid of the occupation variable a by defining the relation 



that, applied in the case of apolar representations, allows the replacement 
of rf by Tj. This new constraint can be interpreted as "if there is a leg 
somewhere in the lattice r , there must also be a particle in the corresponding 
site a. " or, considering water molecules without isolated hydrogen or oxygen 
atoms, a more physical interpretation is "if an hydrogen atom is found, then 
there is at least one oxygen atom attached to it". The hamiltonian of equation 
[T]then becomes 
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With this, the r variables now carry all information about both short- 
range(van der Waals) and long-range (hydrogen bonding) interactions. The 
new expression has a larger number of terms, but a close inspection will 
reveal that such terms are rather simple, with the advantage of having just 
one variable to represent both types of degrees of freedom. 



3. The Finite Lattice Approximation 

In order to obtain an expression for the density of the system, the grand 
partition function was approximated evaluating 



E 
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by enumerating and summing it over all possible configurations |29|, |30 
with (m,n,o,p) denoting the state of the system and E mnop its energy over 
the 2x2 lattice with periodic boundary conditions. 

Subsequently, an expression for the density of the system was obtained, 
and a phase diagram with coexistence curves and a maximum of density was 
generated. 

As it was presented, the partition function was obtained by expanding the 
hamiltonian of equation [Tj on the lattice sites and summing over all states. 
This approach was possible only because the search space in this case was 
severely reduced due to the simplifications applied, allowing the further use 
of a software tool to perform the symbolic summations. The resulting ex- 
pression, 



E 3st = 3 (27 + 108e^ + 72e 2/3{ - v+ ^ + go e 2 ^- 2u+v+ ^ + Ae 3 ^ 2v+ ^ + 8e 3 ^- 4u+2v+ ^ 

+ ]_Q e 4/3(-4M+3u+M) _|_ lQ e 4fS(~3u+3v+fj.) _|_ e 4/3(-2u+3v+/i) _|_ QQ e /3(-Su+6v+3fi) 
+ 3g e /3(-4«+6i>+3^ 

(8) 

combined with the relation 
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Figure 2: Schematics of the four-site lattice used over gray-shaded background. Shaded 
sites outside the gray square correspond to the site inside the square with the same shading, 
indicating the periodic connectivity used. 



1 <91nH _ 1 1 c^S 
gives the normalized expected particle density, 

p3st _ [Q e (2(6u+v+n))P (5 _|_ 4e (4«)/3j _|_ 27 e (16«+M)/3 _|_ 3e (4u+6<;+3/ i )/3 ^ + e ( 4 «)/3) 

(1 + 7 e (4«)/3 + e (8u)/3) + e (4(3t.+M))/3 (io + I6e (4u)/3 + e^)] x 

^ 18e (2(6«+^+M))/3 (5 _|_ 4 e (4«)/3j _|_ 4 e (4«+6^+3^)/3 _|_ e (4«)/3j ^ + 7e (4«)/3 + e (8u)/^ 

+e (4(3«+M))/9 (io + 16e (4u)/3 + + 27e (16u)/3 (l + 4e (M) ' 9 )]~ 1 . 

(10) 

Where \x is the chemical potential, = 1/ksT and V = L 2 is the number 
of sites in the lattice. 

The resulting density for u = 1, v = 1, plotted against chemical potential 
for some values of temperature, depicts curves with the same equilibrium 
densities already found in more complicated models [20l|. as can be seen in 
figure [3j 

The pressure was then obtained by numerically integrating the Gibbs- 
Duhem equation: 

SdT -VdP + Ndfi = (11) 
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and keeping the temperature constant so the first term vanishes. The result- 
ing pressure isotherms can be seen in figure [3J From here we start using the 
overhead bar (e.g. T) to denote dimensionless variables. 




Figure 3: Density isotherms p against chemical potential \x from equation 1101 Triangles 
correspond to T = 1(T 4 , circles to f = 0.1, X's to T = 1.0, and squares to T = 10.0. It is 
possible to notice the presence of three stable phases when T < 0.1. As T gets closer to 
1, both full and empty phases remain, but the kagom phase (p = 0.75) gets unstable and 
vanishes as T increases. If we go further towards high T, we see that p = 1 also doesn't 
last anymore. 



From the pressure curves obtained, it is possible to sketch the phase space 
of the system. Each level of the pressure steps produces a coexistence curve 
along T, in this case, separating regions of gas and low density liquid phases, 
and also between low and high density liquid phases, respectively. 

Investigating the suitable region of the density inf-jU surface, we were 
also able to identify a region of maxima, as can be seen in the resulting phase 
diagram of figure HI 

We notice that, even with the known limitations of finite-lattice approx- 
the results from such small lattice dimensions still allowed 

and 



imations 311 32 
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the occurrence of the same stable phases already shown in 
resulting in a phase diagram very similar to that found by monte carlo sim- 
ulations with larger lattice sizes. 
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Figure 4: Comparison of coexistence curves from Monte Carlo simulations (triangles) [19( 
and finite enumeration (stars) for the transitions Gas-LDL (below), LDL-HDL (above), 
and maximum of density (right). Parameters u and v are set both to 1. 



4. Proposed Density 

Considering the complexity of the expression for the density of the system 
previously obtained by exact calculations (equation [10]), the possibility of 
representing these curves with a simpler expression was considered, so that 
it could be used to phenomenologically approximate the system behavior 
and reproduce some of its features - for example, the anomaly of the density. 
Another important motivation for seeking a simpler solution, is the difficulty 
to perform an analytical integration in the resulting expression of the exact 
density to obtain the pressure, for example. With a simplified and integrable 
representation for the density, it would be possible to obtain expressions 
for physical observables that provide insights about the real behavior of the 
system in study. After some simple calculations, we will detail the theoretical 
hypothesis that motivates this section. 

Removing all the inter-particle interaction terms of the hamiltonian, we 
reduce our model to that of a lattice gas where each particle interacts exclu- 
sively with the adsorbant (the lattice), giving 
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7~Lads — — 




(12) 



So the grand partition function can be easily obtained by exact state 
enumeration in a finite lattice, 



The factor 81 multiplying the partition function in equation [13] arises 
because the number of states per particle powered to the number of sites is 



The following hypothesis was then raised: Since the exact result obtained 
by excluding the short and long-range interaction terms from the model has 
a logistic form, could it be possible to represent the densities of the stable 
liquid phases of the system as components of a superposition of logistic func- 
tions, whose forms are qualitatively compatible with the stable regions of the 
density pi In this representation, what characteristics of the original system 
would be preserved? And what would be the minimum number of terms 
necessary to provide an "acceptable" result from the perspective of structure 
of the phase diagram? 

Thus, an initial, phenomenological expression for the density of the sys- 
tem, consisting of a weighted superposition of partial densities based on the 
equation [TH was obtained, for the particular case where u/v ~ 1. This spe- 
cific case greatly simplifies the modelling, and keeps the region where the 
two liquid phases coexist, because of their importance in the formation of 
the maximum of density that we intend to study. 

Assuming that the anomalous density of our model can be approximated 
by a composition of terms similar to the equation [TH for each phase that 
composes the complete density of the system, the resulting equation is 




(13) 



and the density of the system is 




(14) 



81. 




(15) 
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where p\ are the local variations in the stable partial densities for each 
critical /i c , and //j is the smallest value of the chemical potential where the 
i-th contribution manifests itself, taking into account the existence of the 
densities present in smaller chemical potentials. 

Analyzing the results from simulations and from the exact analytical so- 
lution for the system in the finite lattice with 3 states, we particularize the 
equation [15] for the expression [TB] below, 



+ 



1 _|_ e k B T 1 + e k B T 



(16) 



where p, g i ~ —2 and ftih ~ 2 are the chemical potentials of the gas-LDL 
and LDL-HDL transitions respectively. Such parameteres will depend on the 
fluid model adopted, and can be estimated by simulation. 

The equation [16] is our proposal for the simplified density of a two- 
dimensional liquid polymorphic fluid, in an expansion consisting of the linear 
combination of two logistic functions. We then start working with equation 
TBJfor the normalized particle density in order to measure its ability to ap- 
proximate the behavior of the stable phases of the system. Notice that, in 
this representation, the parameters u and v are implicit. 

In the figure we can see the behavior of the density in equation [TBI as 
a function of the chemical potential /x, for different temperature ranges, in 
comparison with the exact solution of finite lattice for the model of 3 states 
in 4 sites previously shown, in identical conditions of p, and T. 

According to figure [5] the proposed expresssion for density curves gener- 
ates results qualitatively similar to that obtained by direct integration of the 
partition function in the finite lattice. For the case of the specific density 
that we want to approximate, originated from first-principles physical model, 
the parameters /i; had to be estimated by the values of the chemical potential 
where the phase transitions happened in monte carlo simulations, taking the 
low temperature range as a reference. 

Using once again the Gibbs-Duhem relation, we can obtain an analytical 
expression for the pressure 



fix 1 



M+M2 



/ r a mi i r „ , 

p= + -k B T {log \e k B T +e fc srj + 3 log ^1 + e fc s T j J . (17) 

where — ^ is an integration constant obtained from the condition that 
the pressure for zero density must also be zero, fix is the chemical potential 
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Figure 5: Comparison among density curves obtained from the exact solutions in 4 sites 
for three different temperatures (circles— i- 0.1 , squares— > 1.0 and x's— > 10.0), and curves 
generated by logistic expansions (empty triangles—!- 0.1, crosses—!- 1.0 and stars— !• 10.0). 
The match between results is higher at low temperatures, decreasing at f « 1.0 and 
increases again as T rises. 

where the first phase transition , gas — > LDL occurs, and /i2 is the chem- 
ical potential where the second phase transition , LDL — > HDL occurs. 
Isotherms generated from the equation [T7] are presented in figure [6j 

The comparison of the resulting phase diagram with a phase diagram 
obtained by exact integration over a finite lattice can be seen at figure [71 
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Figure 6: Pressure isotherms obtained from equation [17] The sharper curve corresponds 
to f = 1(T 4 , until f = 10° for the smoother one. A close similarity to Monte Carlo results 
from 19] can be noticed. 
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Figure 7: Comparison of coexistence curves for the finite enumeration method (stars) and 
logistic expansion (diamonds) . It is possible to see transitions GAS-LDL (below) and LDL- 
HDL (above). The small curves at right correspond to the maximum of density, captured 
in both approaches. Notice how most properties of the original model are captured by the 
approximation proposed in equations I15[ 1161 and 1171 
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5. Conclusions 



The obtained results led us to conclude that finite integrations on a 2 x 2 
lattice were be able to capture, for a 2-dimensional water-like system, most 
of the qualitative properties of the system observed in previous Monte Carlo 
simulations, such as gas, HDL and LDL phases. Although the temperature 
range of the coexistence lines in our approach is shorter (« 20% less), with 
lower limits, mainly due to the smoothing of transitions because of finite 
size effects, the phases originally observed in Monte Carlo simulation are 
still present. The maximum of density also remains in the low- dimensional 



system, because the quasi-cristalline network in the liquid phase |33[ finds a 
basis in the lattice with minimal resolution but still enough to occur. It is 
surprising that such extreme reduction of dimensionality still allows the phase 
diagram to keep most of its characteristics found in Monte Carlo simulations. 

Regarding the proposed (logistic) representation for the system density, 
most of the features of the phase diagram, obtained previously in results from 
both Monte Carlo simulations and from finite system approximation, remain 
present in the phase diagram obtained from the pressure equation integrated 
from the logistic density. The phase diagram generated using the proposed 
density also has the anomaly found in previous results from simulations and 
finite approximation. 

Analyzing the P — T phase diagram with the pressure obtained by inte- 
grating the logistic density represented in this work by a sum of two more 
fundamental adsorption curves, we believe that the anomaly phenomenon 
is the result of two interpenetrating phases, being the less dense a quasi- 
crystalline phase with vacancies, and the low-density a fluid of non-bonded, 
loose particles. As temperature rises, vacancies in the quasi-crystalline liquid 
are filled with free molecules travelling in the fluid, leading to an increase 
in the density, and as the temperature continues rising, the quasi-crystalline 
structure breaks, making the density finally starting to decrease. 

As a future work, we intend to further investigate these effects both by 
simulations and analytically in larger systems, in two and three dimensions, 
so to get a better understanding about the formation rules underlying the 
maximum of density in a real, three-dimensional fluid. 
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